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This  report  summarizes  the  research  results  obtained  in  the  past  three  years,  and  sup- 

rv 

ported  by  -the-urbove  grant  from  the  U.S.  Army  Research  Office.  Only  the  main  results  are 
summarized  and  for  the  details  the  reader  is  referred  to  the  appropriate  publications.  All  refer¬ 
ences  here  refer  to  our  publications,  a  list  of  which  is  included  at  the  end  of  this  report.  Also 
included  at  the  end  of  the  report  are  the  abstracts  of  the  publications. 

The  main  focus  of  this  project  was  the  evaluation  of  high  resolution  spectrum  estimation 
methods  with  particular  emphasis  on  state  space  based  methods.  Of  particular  interest  are 
finite  precision  etrects  (numerical  issues),  and  the  statistical  efficiency  of  high  resolution 
methods.  We  believe  our  results  have  in  a  significant  manner  helped  in  the  understanding  of 
these  methods.  The  results  are  divided  into  and  discussed  under  three  categories:  Spectrum 
Estimation,  Array  processing  and  Floating  Point  arithmetic.  ;  \  : 

,  i 

Spectrum  Estimation: 

In  this  area,  the  problem  of  frequency  estimation  of  sinusoids  from  time  series  is  considered. 
Of  particular  interest  is  the  estimation  of  frequencies  from  short  and  noisy  data  records. 
Fourier  based  methods  have  many  deficiencies  and  are  not  appropriate  in  this  context.  On  the 
other  hand,  model  based  methods  have  proven  to  have  high  resolution  capabilty  and  have 
attracted  a  great  deal  of  interest.  Numei  ■us  model  based  methods  have  been  developed  and 
applied  to  the  problem  of  frequency  estimation.  Of  particular  interest  are  Singular  Value 
Decomposition  (SVD)  based  methods  which  show  great  promise.  Such  methods  are  analyzed 
in  our  work  with  particular  attention  being  paid  to  state  space  based  methods,  i.e.  methods 
developed  using  a  state  space  representation  for  the  models.  The  main  results  in  this  context 
are  as  follows: 

1.  Evaluation  of  model  based  methods  is  fairly  involved,  and  in  a  great  deal  of  the  existing 
work,  computer  simulations  have  been  used  to  support  the  methods.  In  this  research,  an 
insightful  approach  to  evaluate  the  methods  has  been  developed.  It  was  recognized  that  most 
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model  based  methods  can  be  divided  into  two  distinct  steps.  Step  1  consists  of  estimating  the 
parameters  of  the  model  from  the  given  data,  and  Step  II  consists  of  extracting  the  required 
information,  frequency  estimates  in  this  case,  from  the  estimated  parameters.  Both  steps  are 
important  to  the  overall  quality  of  the  estimate.  Step  I  requires  reliable  estimation  of  the 
parameter  set,  and  Step  II  requires  that  the  parameter  set  have  low  parameter  sensitivity.  An 
evaluation  of  both  steps  are  conducted  to  obtain  insight  into  SVD  based  methods.  A  more 
extensive  discussion  of  the  pros  and  cons  of  such  an  evaluation  can  be  found  in  [7,13,2,4,1]. 

2.  The  two  step  analysis  was  first  applied  to  the  SVD  based  method  of  Tufts  and  Kumaresan. 
In  this  method,  the  linear  predictability  of  sinusoids  is  exploited  to  obtain  the  frequency  esti¬ 
mates.  The  procedure  consists  of  first  estimating  the  linear  prediction  parameters  using  a  SVD 
based  approach,  and  then  estimating  the  frequencies  from  the  roots  of  the  linear  prediction 
polynomial.  It  was  shown  that  for  step  II,  the  parameter  sensitivity,  a  large  polynomial  along 
with  the  minimum  norm  solution  are  the  key.  This  allows  the  redundant  roots  of  the  polyno¬ 
mial  to  play  an  important  role  in  reducing  the  sensitivity  of  the  desired  roots  to  perturbations 
in  the  coefficients  [7,2,4].  In  Step  I,  for  reliable  estimation,  it  is  shown  that  the  continuity  of 
the  generalized  inverse  and  the  concept  of  angle  between  subspaces  play  an  important  role. 
The  continuity  concept  helps  explain  the  need  for  a  low  rank  approximation  in  the  estimation 
process,  and  the  quality  of  the  approximation  is  appraised  using  the  concept  of  angle  between 
subspaces  [7,2,4]. 

3.  Step  I,  reliable  estimation  is  further  explored  in  [1]  and  [3,6],  In  [1],  it  is  shown  that  the 
least  squares  problem  solved  in  the  Kumaresan-Prouy  approach  can  be  ill-conditioned  leading 
to  degradation  in  the  quality  of  the  estimates.  In  particular,  this  provided  a  theoretical  explana¬ 
tion  to  a  strange  behavior  that  was  observed  by  others  when  the  method  was  used  for  direc¬ 
tion  finding,  namely  that  for  a  particular  two  source  scenario  the  performance  degraded  as  the 
separation  increased  as  opposed  to  improving  as  is  usually  the  case.  In  [6,3],  the  concept  of 
angle  between  subspaces  is  further  explored.  This  concept  provides  a  procedure  to  decide 
when  SVD  based  methods  fail,  i.e.  gives  an  indication  of  the  theoretical  limits  of  SVD  based 
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methods.  Also  it  was  found  that  existing  SVD  based  model  based  methods  did  not  exploit  the 
information  in  the  singular  vectors  maximally.  In  [6],  a  non-linear  method  starting  from  the 
singular  vectors  was  developed.  The  method  consisted  of  fitting  a  structured  subspace  to  the 
space  spanned  by  the  principal  singular  vectors.  The  subspace  fitting  method  turns  out  to  be 
better  than  the  existing  SVD  methods  but  is  computationally  more  intense.  However,  the  ideas 
developed  provide  important  insight  namely  that  the  existing  methods  do  not  necessarily  make 
the  best  use  of  the  principal  singular  vectors. 

4.  The  state  space  methods  are  studied  in  [2,13,4],  The  approach  differs  from  the  Tufts  and 
Kumaresan  approach  in  the  fact  that  it  is  based  on  a  state  space  parameterization  of  the  model 
instead  of  the  polynomial  parameterization.  The  state  space  methods  use  a  SVD  Kas?d 
approach  to  estimate  a  state  transition  matrix,  and  the  frequency  estimates  are  obtained  from 
the  eigenvalues  of  the  matrix.  The  issues  related  to  estimation,  step  I,  turn  out  to  be  the  same 
as  before  and  so  only  parameter  sensitivity,  Step  II.  is  examined  in  detail.  It  is  shown  in  the 
context  of  sinusoid  frequency  estimation  that  from  a  parameter  sensitivity  point  of  view  the 
best  state  transition  matrix  to  estimate  is  an  unitary  matrix.  Procedures  to  estimate  such  uni¬ 
tary  matrices  are  then  explored.  It  is  shown  that  the  Toepiuz  Approximation  Method  (TAM) 
and  the  Direct  Data  Approximation  (DDA),  two  existing  methods  do  estimate  such  robust 
matrices  in  a  reliable  manner  [2,13.4], 

5.  In  [5],  attempts  to  further  improve  the  estimation  procedure  are  made.  The  fact  that  the 
matrix  obtained  using  TAM  is  ideally  unitary  is  incorporated.  This  leads  to  a  constrained  least 
squares  problem  with  a  fairly  simple  solution  [5], 

6.  Interesting  connections  between  state  space  methods  and  the  Matrix  Pencil  approach  are 
shown  in  [18],  It  is  shown  that  the  state  space  approach  and  the  Matrix  Pencil  approach  are 
very  closely  related  [18]  in  that  if  the  state  space  method  can  be  viewed  as  estimating  a 
matrix  F=AB ,  then  the  matrix  pencil  method  amounts  to  estimating  F  =BA .  Note  that  the 
nonzero  eigenvalues  of  AB  and  BA  are  the  same  leading  to  the  same  frequency  estimates.  This 
implies  that  the  performance  of  the  two  methods  are  nearly  identical.  The  superiority  of  State 
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Space  methods  over  the  polynomial  approach,  e.g.  Tufts  and  Kumaresan  method,  follow  from 
the  results  of  Hua  and  Sarkar  who  show  that  Matrix  Pencil  methods  are  usually  superior  or 
comparable  in  performance  statistically  to  polynomial  methods.  In  addition  to  the  statistical 
performance,  state  space  methods  are  often  computationally  more  desirable.  One  particularly 
advantageous  feature  of  state  space  methods  is  that  extracting  the  frequency  estimate  from  the 
model  parameters  is  much  simpler. 

7.  In  addition  to  the  sinusoid  problem,  the  state  space  formulation  has  wider  applicability  and 
these  are  highlighted  in  [2,20].  The  applications  of  state  space  ideas  to  variety  of  problems, 
e.g.  estimating  parameters  of  damped  exponentials,  ARMA  power  spectrum  etc.  are  discussed. 
A  key  outcome  of  this  research  is  evidence  to  indicate  that  the  state  space  formulism  offers 
numerous  advantages.  For  instance,  it  provides  a  convenient  framework  for  exposing  struc¬ 
ture  inherent  in  signal  processing  problems.  Also  the  flexibity  it  provides  in  the  parameter 
domain  can  be  exploited  to  lower  parameter  sensitivity.  Furthermore  reliable  methods  exist  to 
estimate  these  robust  parameter  sets.  All  these  attributes  are  very  beneficial  to  signal  process¬ 
ing  algorithms. 

Snbspacc  Based  Array  Processing: 

Similar  to  methods  in  spectrum  estimation,  methods  based  on  SVD  or  Eigen  decomposition 
are  popular  for  the  problem  of  direction  finding  or  the  direction  of  arrival  (DOA)  estimation 
problem  using  a  sensor  array.  In  this  problem  the  data  consists  of  the  array  output  at  different 
time  instants.  The  output  of  the  array  at  any  time  instant  is  referred  to  as  a  snapshot.  Using 
the  snapshots,  an  estimate  of  the  covariance  matrix  of  the  data  is  obtained.  Then  an  estimate 
of  the  DOA  is  obtained  by  using  an  eigendecomposition  of  the  covariance  matrix.  Such 
methods  are  often  referred  to  as  Subspace  based  methods.  There  are  numerous  subspace  based 
methods  and  the  purpose  of  this  research  was  to  evaluate  their  pros  and  cons.  Our  research 
analyzed  some  of  the  popular  methods,  namely  MUSIC,  ESPRIT,  TAM  and  the  Minimum- 
Norm  method,  under  the  assumption  that  the  snapshots  were  independent  and  that  the  data 
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was  circularly  Gaussian.  Also  the  estimate  of  the  covariance  matrix  considered  in  the  analysis 
is  that  obtained  by  taking  the  average  of  the  outer  product  of  the  snapshot  vectors.  The  main 
results  of  the  analysis  are  as  follows: 

1.  It  is  shown  that  for  the  linear  equispaced  sensor  array  case,  ESPRIT  can  be  treated  via  a 
state  space  formulation.  The  matrix  obtained  using  ESPRIT  is  shown  to  be  related  to  the 
matrix  obtained  using  TAM,  by  a  diagonal  transformation.  TAM  and  (Least  Squares) 
ESPRIT  are  shown  to  have  identical  statistical  properties.  Also  Total  Least  Square  ESPRIT 
and  Least  Squares  ESPRIT  are  shown  to  have  the  same  asymptotic  mean  squared  error. 
Closed  form  expressions  are  derived  for  the  mean  squared  error.  Simple  closed  form  expres¬ 
sions  are  derived  for  the  one  and  two  source  case  leading  to  interesting  insights  [8,14]. 

2.  For  the  linear  equispaced  array  case.  Closed  form  expressions  for  the  asymptotic  mean 
squared  error  in  the  estimate  of  the  direction  of  arrival  for  the  Minimum-Norm  method  and 
Root-Music  are  also  obtained  [12,17,9,15],  Root  Music  is  a  variant  of  MUSIC  in  that  the 
DOA  estimates  are  obtained  from  the  roots  of  a  polynomial  obtained  by  performing  a  spectral 
factorization  of  the  null  spectrum  [9,15].  It  is  shown  that  the  DOA  estimates  obtained  using 
Root-Music  has  a  smaller  mean  squared  error  compared  to  the  Minimum-Norm  method. 
Also  the  estimates  obtained  using  Root-MUSIC  become  significantly  better  as  the  array  length 
increases.  Simple  closed  form  expressions  are  derived  for  the  one  and  two  source  case.  Simi¬ 
lar  conclusion  also  apply  when  Root-MUSIC  and  ESPRIT  are  compared.  Root-MUSIC 
appears  to  be  the  best  of  the  three  under  the  assumed  conditions. 

3.  It  was  shown  that  in  the  case  of  Root-MUSIC,  though  the  errors  in  the  DOA  estimates 
were  small,  the  error  in  the  roots  was  large.  This  implies  that  for  Root-Music  the  errors  in 
the  roots  were  mainly  radial.  This  observation  is  fairly  important  in  that  it  implies  trouble 
when  using  a  spectral  approach.  The  radial  errors  result  ir  less  well  defined  peaks.  If  the  roots 
are  close  together,  due  to  the  large  radial  error  there  may  not  even  be  two  peaks  in  the  spec¬ 
trum.  This  in  turn  results  in  a  apparent  loss  in  resolution  and  a  higher  threshold  SNR  for 
(spectral)  MUSIC  (compared  to  the  spectral  Minimum-Norm  method).  Also  it  is  shown  that 
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asymptotically  there  is  no  difference  between  using  a  spectral  approach  versus  using  a  rooting 
procedure.  However,  for  the  reasons  explained  above,  rooting  is  desirable  for  MUSIC  wher¬ 
ever  possible. 

In  summary',  the  above  results  provide  valuable  insight  into  the  above  methods  in  the 
context  of  a  linear  equispaced  sensor  array  and  when  the  sources  are  uncorrelated  or  partially 
correlated.  To  overcome  the  problems  that  arise  when  the  sources  are  fully  correlated,  spatial 
smoothing  is  employed.  In  this  approach  the  array  is  divided  into  smaller  (overlapping)  subar¬ 
rays  and  the  overall  covariance  estimate  obtained  by  averaging  over  the  covariance  estimate 
obtained  from  each  su barmy  [19].  The  main  results  in  the  context  of  spatial  smoothing  are  as 
follows. 

4.  The  statistics  of  the  eigenvectors  of  a  spatially  smoothed  forward-backward  covariance 
matrix  are  presented  [19].  This  result  is  of  general  interest,  and  is  used  to  evaluate  subspace 
based  array  processing  methods.  The  procedure  used  to  obtain  these  statistics  consists  of  using 
a  first  order  expansion  of  the  signal  space  eigenvectors.  The  approach  is  fairly  general  and  can 
be  easily  generalized  to  deal  with  other  scenarios,  i.e.  effect  of  sensor  perturbations,  deter¬ 
ministic  data  model  etc  [19]. 

5.  The  statistics  are  used  to  study  the  effect  of  spatial  smoothing  on  Root-MUSIC  and  the 
Nlii.iiuum-Nor.r.  ...vthod.  Expressions  for  the  mean  squared  e^r  in  the  DOA  estimates  are 
derived  [19],  The  methods  are  compared  and  it  is  found  that  spatial  smoothing  is  more 
beneficial  to  the  Minimum-Norm  method  compared  to  MUSIC.  In  our  earlier  work  it  was 
shown  that  when  no  smoothing  was  used,  MUSIC  was  mper-inr  m  the  Minimum-Norm 
method  [12,17],  Here  it  is  found  that  by  properly  choosing  the  number  of  subarrays,  the  per¬ 
formance  of  the  Minimum-Norm  method  can  be  made  comparable  to  MUSIC  [19].  In  the 
case  of  MUSIC,  though  the  DOA  estimates  usually  tend  to  deteriorate  as  the  number  of  subar- 
ravs  increase,  the  mean  square  error  in  the  roots  improve.  Hence,  if  spectral  MUSIC  is  used, 
the  peaks  will  be  sharper  as  the  number  of  subarrays  is  increased  giving  the  false  impression 


of  better  resolution. 
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6.  The  effect  of  spatial  smoothing  on  ESPRIT  is  also  studied.  Similar  to  the  Minimum  Norm 
method,  spatial  smoothing  is  beneficial  to  ESPRIT  and  significantly  improves  the  estimates. 

7.  Some  implementation  (computer  architecture)  issues  are  also  considered  in  [10].  It  is  shown 
that  since  State  Space  based  DOA  methods,  and  ESPRIT  essentially  require  only  matrix 
operations,  they  are  very  suitable  for  systolic/wavefront  implementation  [10,20]. 

Finite  Precision  Floating  Point  Arithmetic: 

The  effect  of  finite  precision  floating  point  arithmetic  in  digital  signal  processing  is  not 
well  understood.  However  with  the  growing  availability  of  floating  point  chips,  understanding 
the  effect  of  finite  precision  floating  point  arithmetic  is  an  increasingly  important  problem. 
Our  work  in  this  area  is  very  fundamental,  and  has  the  potential  of  providing  insight  into 
many  applications,  e.g.  effect  of  finite  precision  floating  point  arithmetic  on  digital  filters 
[11,16],  effect  of  finite  precision  floating  point  arithmetic  on  adaptive  algorithms  etc.  The 
main  results  are  as  follow's: 

1.  As  a  starting  point  for  the  study,  the  problem  of  digital  filters  is  addressed  with  particular 
emphasis  on  robust  state  space  digital  filters.  In  this  context,  fixed  point  arithmetic  has  been 
extensively  studied.  The  limited  success  with  analyzing  the  effect  of  floating  point  arithmetic 
on  digital  filters  has  been  due  to  the  fact  that  the  roundoff  errors  are  correlated  with  the  sig¬ 
nals  making  the  analysis  tractable.  In  our  work  to  overcome  this  difficulty,  a  very  unique 
interpretation  of  floating  point  arithmetic  is  provided.  It  is  shown  that  the  inner  product  a'x. 
an  operation  basic  to  linear  time  invariant  systems,  when  computed  using  finite  precision 
floating  point  arithmetic  is  equivalent  to  (a-Ai)'x,  t.e.  the  exact  inner  product  of  a  perturbed  a 
and  .r.  This  interpretation  in  the  context  of  digital  filters  implies  that  the  roundoff  noise  can  be 
determined  by  perturbing  the  filter  parameters.  It  provides  a  strong  link  between  coefficient 
sensitivity  and  roundoff  noise  and  also  makes  the  mathematics  much  more  tractable. 

2.  An  expression  for  the  variance  of  the  output  roundoff  noise  is  derived  and  shown  to  depend 
on  the  filter  parameters,  in  particular  the  Observabilty  and  Controllabiltv  Grammians.  Unlike 


-  9  - 


in  the  fixed  point  ease,  the  variance  also  depends  on  the  input  signal  statistics.  For  the  case 
where  double  pree: .  jn  accumula'ion  is  used,  it  is  shown  that  the  optimal  filters  are  similar  to 
those  obtain  '  .vhen  using  fixed  point  arithmetic. 
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ABSTRACT 

In  t his  paoer.  we  perform  a  numerical  analysis  of  the 
Kumaresan-Pronv  method,  a  method  for  estimating  the  frequencies  of 
sinusoids  m  white  noise.  Like  most  model  based  methods,  in  this 
soproach  the  re'iafeiJitv  of  the  Estimation  procedure  and  Parameter 
sensicivicv*  are  central  :o  determining  ns  performance.  We  6m  exam¬ 
ine  the  pnorameter  «ensmvicv  issue  and  explain  how  the  excess  filter 
length  helos  combat  numertc.ai  sensitivity  We  examine  the  estimation 
procedure  and  mnic.ate  how  *he  conceot  of  angle  between  subsoaces 
can  be  used  'o  determine  ’.he  quality  of  estimation.  An  interesting 
outcome  of  'he  analysis  is  the  explanation  provided  to  the  limitation 
of  the  Kumar-san-Prony  algorithm  when  applied  to  the  direction 
finding  problem. 


L  INTRODUCTION 

In  this  paoer  we  perform  a  numerical  analysis  of  the 
Kumarcsan.pronv  (KP1  method.  The  KP  method  is  a  special  case  of 
the  general  Tufts  and  Kumaresan  appraoch.  a  SVD  based  linear  pred¬ 
iction  method  for  *snmacing  the  fr*quencies  of  sinusoids  in  noise  tl|. 
To  begin  the  analysts  we  note  that  the  model  based  methods  can  tn 
general  be  divided  into  two  distinct  steps  as  in  |2J.  Step  I  con¬ 

sists  of  estimating  i  parameter  set  that  describes  the  model  And  step 
II  consists  of  determining  the  relevant  information  from  the  parame¬ 
ter  -scimates.  Both  steps  are  important  to  the  overall  mccess  of  the 
metnod.  lll-condicionmg  \t  either  step  can  adversely  effect  the  overall 
performance  of  the  method  and  needs  to  be  avoided  12.21. 

For  the  sinusoid  problem  let  us  denote  the  parameter  set  bv 
9  '■  '  &\  .  ’  ind  the  quantity  of  interest  are  the  frequencies 

w,  f.  An  -rror  A 6  in  che  '•stimate  of  9  results  in  a  error  Aw,  in  the 
frequency,  where 


-  E  —  Ay.  ■=*  Awr  AW 

where 

Aw,  ,  ck* /,  | 

d0,  '  59,  off,  j 

and  A0  ’.s  denned  similarly 


S  represents  the  sensitivity  of  the  parameter  set  and  /?,  the  estima¬ 
tion  error  covariance.  In  step  I,  one  attempts  to  minimue  ;|  R, ,  |+ 
However  for  the  overall  success  of  the  method.  ill-conditioning  in  step 
II  has  to  be  avoided,  i.e.  5  has  to  be  small.  Reduced  sensitivity  also 
reduces  the  demands  on  the  estimation  process.  Though  an  optimal¬ 
ity  in  both  steps  is  desirable,  in  practice,  this  mar  be  hard  to  achieve. 
However  it  is  clear  that  an  analysis  of  the  two  steps  separately  can  be 
useful  in  understanding  and  reaching  the  best  compromise.  In  this 
paper  we  conduct  such  %  two  step  analysis  of  the  KP  method. 


n.  BACKGROUND 

The  observations  y (k)  consists  of  the  sinusoidal  <ignai  r{t)  cor¬ 
rupted  by  white  noise  a(ifc),  t.e. 

y{k)  -  z[k)  *■  rt(*\  ,  k  -  l  jU; 

The  signal 


r{k) 


{ lb> 


where  c,  ’ t  are  the  complex  amplitudes  and  w,  ’.he  fluencies. 

A  popular  procedure  to  estimate  the  frequencies  is  to  use  an 
Linear  Prediction  (LP)  approach  4i.  In  the  absence  of  noise  the 
sinusoidal  signaj  is  exactly  predictable,  i.e. 


(2) 


The  teros  of  the  polynomial  formed  from  the  coefficients 

p 

ff(*)  - 1  -  £**-■  .  gives  an  estimate  of  the  frequencies.  Using  the 

»-i 

LP  approach  the  problem  ol  frequency  estimation  is  reduced  to  that 
of  estimating  the  coefficients  9,.  For  the  estimation  of  the  coefficients 
«  least  square  approach  can  be  used.  In  the  presence  of  noise  the  pro¬ 
cedure  is  less  efficient  and  improvements  can  be- achieved  by  nsin*  a 
forward  backward  prediction  approach  i5,6l.  It  can  be  mown  that  the 
same  coefficients  appear  in  the  backward  prediction  equation,  i.e. 

y(l)  -  £ ;.’y(*+0  -  0. 


£  '  Aw.i  ’  -  Aw'  5|a0  A0'j  Aw  <-  5  ||  R, J| 
where 


5 


Aui  f 


V 


W 


and  Rt 
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"l.'nle.s  mention'd  otnerwi.e  ihe  norm  referred  to  in  this  paper  is 
'.he  ,  norm 


Using  the  forward  backward  approach,  the  problem  reduces  to  solving 
the  following  set  ol  linear  equal  ions. 


'j(p) 

y|p-i) 

y|D 

yip  -l) 

yip  - 1) 

y(p  -2| 

y\-l 

y(p-i| 

9\ 

yiiV-l) 

y('V~2| 

vi'V-n) 

yl.VI 

ir'l'2) 

y  12) 

-  y  i  p  -  1 1 

* 

y  1 1) 

y'(3) 

y'f4) 

y’(p*N 

y'lN 

9* 

rl.v-iM-i) 

V  *f  ,V-p  wn) 

yl.VI 

y  ’ !  ,v  -  o ) 
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ABSTRACT 

In  this  paper.  •-*  e  examine  the  issue  of  parameter  sensitivity  in  spec¬ 
trum  estimation.  !n  rarncular.  state  space  models  are  considered  and 
robust  coordinate  <\<irms  for  «pertr«rn  estimation  are  identified.  F^r 
the  sinusoid  problem  a  is  shown  that  the  ideal  parameter  set  involves 
eciim/iun?  i  umtarv  matrix  and  simple  procedures  to  estimate  such  a 
matrix  are  identified.  T  »r  the  damped  sinusoid  and  the  general 
ARM  A  apectrum  estimation  problem  it  is  shown  that  the  baianced 
coordinate*  are  robust  and  reliable 


I.  INTRODUCTION 


5  represents  the  sensitivity  of  the  parameter  «et  and  !it  he  estima¬ 
tion  error  covariance.  In  step  1,  one  attempts' to  minimize 
However  for  (he  overall  success  of  »he  method,  ill-conditioning  n  «tep 
II  has  io  be  avoided,  i  e.  .?  has  in  fie  small  Reduced  sensitivity  also 
reduces  the  demands  on  the  estimation  process.  Though  an  optimal¬ 
ity  in  both  steps  is  desirable,  in  practice,  this  may  be  hard  achieve 
fn  this  paper  we  consider  the  prooiern  f  sinusoids,  damped 
sinusoids  and  (hen  (he  general  rational  modeling  (ARMA)  problem. 
In  all  these  cases,  we  concentrate  on  <tep  II.  the  parameterization 
issue  and  attempt  to  determine  robust  parameter  sets.  For  this  pur¬ 
pose  we  focus  our  attention  on  <taie  space  parameterization  and  pay 
particular  attention  to  balanced  coordinate  systems  !4,5.6j. 


In  recent  vrars  therp  Has  h*en  a  great  deal  of  interest  in  mode! 
based  spectrum  estimation  methods  ill.  They  provide  a  tool  to  extend 
i  he  data  beyond  the  observation  interval  and  are  capable  of  providing 
high  resolution  spectral  estimates.  In  this  paper  we  conduct  a  numeri¬ 
cal  stud  v  of  -.ome  spectral  estimation  methods  with  particular 
emphasis  on  slate  space  modeling.  To  begin  the  analysis  we  note  that 
the  model  based  methods  can  in  general  be  divided  into  two  distinct 
steps  as  in  fig  (1)  -21  Step  I  consists  of  estimating  a  parameter  set 
that  describes  the  model  and  step  II  .-onsists  of  determining  the 
relevant  information  from  the  parameter  estimates.  Both  steps  are 
important  to  the  overall  success  of  the  method.  Ill-conditioning  at 
either  step  ran  adversely  ‘•tfert  the  overall  performance  of  the  method 
and  neHs  »o  be  avoided  '2.3l. 

For  ”\ampie  in  the  sinusoid  problem,  let  us  denote  the  parame- 
*  er  set  by  B  -  •  ^  .  u, ,  .9,  I ,  and  the  quantity  of  interest  are  the 

frequencies  w,  i.  An  error  A0  in  the  estimate  of  B  results  in  \  error 
A  a/,-  in  the  frequency,  where 


Ad.  =  A-/  AO 


an<i 


A^ 


Au/f 


3 w, 


with  A(9  being  :»*fm"d  similarly  Here  we  will  use  ’’  £  '*  to  denote 
transpose,  '*  *  "  to  denote  complex  conjugate,  and  '*  //  "  to  denote 
•ornpirx  conjugate  -  ran«ft.«» 

K  A-..  ’  A-  .  1  AB  Al-V'l  A*.  .1 


A-  . 


■V  ,  I 

—  ■  in.|  |A0  AB" 


'This  research  supported  in  part  by  the  ARMY  Research 

Office  under  Grant  No.  D  AAL-03-86-K-0107  and  the  the 
National  Science  Foundation  under  Grant  No.  ECS-^d03732 
U."nie«x  meruione't  otherwise  >he  norm  used  is  the  ’-norm. 


II.  SINUSOID  PROBLEM 

For  'he  sinusoid  problem,  the  data  is  the  sum  of  complex 
exponentials,  i.e. 

y(*l  -  £  r,  t ' 

.  -  I 

Such  a  signal  can  be  modeled  by  a  state  space  model 

r...  -  F  v,  (i») 

y (*)  -  h  Y,  .  (Ib) 

where  the  eigenvalues  of  F  are  of  unit  amplitude  and  equal  to  t'  '  J7). 
One  interesting  realization  is  when 


F  =*  dxaq  1 1 


(2) 


and  r  (0)*  *  (c,  ,Cj,  .  ,  cp )  and  h  «  ( 1.1 ,...,  1).  The  triple  f  F .z  (0),h ) 

characterize  the  model  and  are  not  unique.  If  (F*.r{0).A|  is  one 
representation  then  (  7*“*  FT.  T~[  :  (0),h  T]  is  also  a  represer*  ation  for 
any  nonsingular  matrix  T.  Here  since  the  eigenvalues  of  F  are  of 
interest,  a  matrix  whose  eigenvalues  are  insensitive  to  perturbation  is 
desirable  for  robustness  in  step  II.  This  issue  is  examined  in  detail 
later.  Now  we  present  some  factorizations  that  can  be  used  to  esti¬ 
mate  the  state  space  parameters. 

From  the  state  space  equations,  for  the  sinusoidal  signal  it  can 
be  shown  rhat 


y(*  )  •*  hfk  z  (0). 

The  data  Hankel  matrix  can  be  factorized  into 


(3) 


y(D) 

y(U 

y(’l 

h 

yC) 

y  1 2 ) 

y(’) 

hr 

y(-l 

y(’l 

vCI 

hF: 
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ABSTRACT 

In  this  paper,  we  examine  the  problem  of  evaluating  the 
quality  of  low  rank  approximants  obtained  using  the  Singular 
Value  Decomposition  (SVD).  a  commonly  used  procedure  in 
the  recently  popular  Subspace  Based  methods.  The  concept  of 
angle  between  subspaces  is  used  to  quantify  the  nature  of  the 
approximant.  This  measure  allows  one  to  determine  when  a 
SVD  based  procedures  fails,  and  also  suggests  procedure  to 
maximize  the  perrormance  of  tnese  methods.  The  potential  use 
of  this  result  is  demonstrated  by  applying  it  to  the  problem  of 
esnmaung  the  frequencies  of  mulnpie  sinusoids  in  noise. 


I.  INTRODUCTION 

Singular  Value  Decomposition  (SVD)  has  become  a  very 
popular  tool  in  signal  processing  problems.  For  instance  it  is 
used  widely  in  Least  Squares  problems  which  arise  naturally 
in  many  signal  processing  problems.  In  such  problems  one 
seeks  a  vector  g  which  best  solves 

'  I  =h. 

where  A  is  a  mxn  matrix,  i  is  a  nx\  vector  and  A  is  a  mxl  vec¬ 
tor.  Of  particular  interest  is  the  case  when  the  rank  of  matrix  A 
is  p  where  p  <  minon./i).  The  best  minimum  norm  soluuon  to 
this  least  squares  problem  is  given  by 

p  =  A  *  h  . 

where  a'  denotes  the  pseudoinverse  of  A  and  is  often  com¬ 
puted  using  a  SVD  (Singular  Value  Decomposition)  [1.2). 
Normallv  A  and  A  are  inexact  due  to  the  presence  of  noise.  Let 
the  perturbed  matrix  and  vector  be  denoted  by  s  and  b  respec¬ 
tively.  i.e. 

B  -  A  E .  and  b  -  h  +  e 

Usually  matrix  B  will  be  full  rank,  i.e.  rank  of  8  is  equal  to  t 
where  r  -  minim  zi).  Under  these  conditions  an  often  used  pro¬ 
cedure  is  to  replace  8  by  a  p  rank  approximant  B  before  taking 
the  pseudo  inverse.  An  effective  tool  for  such  an  approxima¬ 
tion  ts  provided  by  the  Singular  Value  Decomposition  (SVD) 
[  1,2).  If  the  SVD  of  B  is  given  by 


where  V  and  f  are  the  left  and  right  singular  vectors,  and  I  is  a 
diagonal  matrix  containing  the  singular  values  in  descending 
order,  i.e.  rr,.a7.  .o,.  The  approximant  is  obtained  by  retain¬ 


ing  the  first  p  principal  singular  vectors  (  l/,andi';  )  and  princi¬ 
pal  singular  values  (I,),  i.e. 

B  =U,Z,Vf. 

and 

=  rr1  UHU 

The  least  square  solution,  denoted  by  x.  is  then  determined  as 
x  =  B*b  . 

There  are  a  number  of  examples  in  signal  processing 
which  can  be  treated  in  the  above  framework,  e.g.  estimating 
the  frequencies  of  multiple  sinusoids  [3,4],  estimating  the 
direction  of  arrival  [5.6|  etc.  SVD  is  also  used  extensively  in 
approximating  a  given  matrix  by  a  low  rank  matrix  and  in 
some  application  the  subspaces  associated  with  the  singular 
vectors  are  used  for  parameter  estimation  [12.9).  Such  an 
application  will  be  pursued  in  sec  III. 

Inspite  of  the  success  of  SVD  based  procedures,  the  ques¬ 
tion  of  evaluating  the  estimation  procedure  still  remains.  Some 
important  questions  that  still  needs  to  be  addressed  are  the  fol¬ 
lowing: 

a)  How  good  is  the  approximant  obtained  by  using  the  SVD  for 
estimation  problems? 

b)  When  do  such  SVD  based  approximation  procedures  fail  to 
be  advantageous  or  useful. 

These  are  the  issues  that  we  attempt  to  address  in  this  paper. 
Some  results  were  first  presented  in  [7.8.9). 


II.  MAIN  RESULT 

For  evaluating  the  quality  of  the  approximant.  and  hence 
the  reliability  of  the  estimation  procedure,  the  concept  of  angle 
between  subspaces  is  considered  in  detail  here.  The  angle 
between  two  subspaces  L  and  M  is  defined  as 

nSm9(LM'ilU=  ll(/  -Pu)  Pi  lb=  »ri  PL  lb 

where  P*  and  PL  are  the  projection  operators  onto  the  sub¬ 
spaces  M  and  L  respectively,  and  is  a  natural  generalization  of 
the  angle  between  vectors  (2.10,11).  If  the  dimension  of  the 
subspaces  are  equal,  then 

i ir* -PL  n  =  ii ri  rL  ii  =  ii rtru  n. 

Consequently,  the  concept  of  angle  between  subspaces  pro¬ 
vides  us  a  useful  tool  to  characterize  the  difference  between 
two  subspaces.  It  is  desirable  to  have  ll P«  -  rL  tl  <  I.  since  this 
ensures  that  no  vector  in  the  subspace  M  is  orthogonal  to  a  vec¬ 
tor  tn  the  subspace  L  [2.10.11). 
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ABSTRACT 

Ln  this  paper,  we  consider  the  problem  of  determining  reliable  model  based  procedures 
for  estimating  the  frequencies  of  sinusoids  from  noisy  measurements.  This  is  done  by  analyz¬ 
ing  the  two  steps  involved  in  such  a  model  based  approach.  One  of  the  steps  involves  analyz¬ 
ing  the  parameter  sensitivity  of  the  model  parameter  set.  and  the  other  involves  examining  the 
reliability  of  the  estimation  procedure  used  to  determine  the  parameter  set  For  this  study,  we 
examine  the  use  of  state  space  models  for  estimating  the  frequencies  of  multiple  sinusoids  in 
noise.  In  particular,  robust  parameter  sets  are  identified,  and  procedures  to  estimate  the  set  are 
presented.  It  is  shown  that  the  ideal  parameter  set  involves  estimating  an  unitary  matrix,  and 
then  computing  its  eigenvalues  to  obtain  an  estimate  of  the  frequencies.  Reliable  Singular 
Value  Decomposition  (SVD)  based  procedures  to  estimate  unitary  or  near  unitary  matrices 
from  covariance  and  tune  series  data  are  presented. 


L  INTRODUCTION 

The  ability  of  model  based  spectral  estimation  procedures  to  produce  high  resolution  esti¬ 
mates  from  finite  data  records  has  attracted  a  great  deal  of  attention  [1-3].  A  multitude  of  tech¬ 
niques  have  been  suggested  and  for  finite  data  records  comparisons  based  on  simulations  have 
been  dene.  Ln  this  paper,  a  two  step  procedure  is  uccd  to  study  the  numerical  issues  involved 
in  such  model  based  methods.  The  main  motivation  for  such  an  approach  is  the  fact  that  model 
based  methods  can  in  general  be  divided  into  two  distinct  steps  as  in  fig.(V).  Step  1  consists  of 
estimating  a  parameter  set  that  describes  the  model,  and  step  II  consists  of  determining  the 
relevant  information  from  the  parameter  estimates.  Both  steps  are  important  to  the  overall  suc¬ 
cess  of  he  method.  Hi-conditioning  at  either  step  can  adversely  effect  the  overall  performance 
of  he  method  and  needs  to  be  avoided.  Reliable  estimation  (Step  I]  is  dearly  an  important 
issue  and  has  been  given  a  fair  amount  of  attention.  For  example,  in  he  dctcrmmisuc  least 
squares  problem  reliable  estimation  procedures,  c.g.  QR  approach.  SVD  approach,  have  been 
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ABSTRACT 

In  this  paper,  we  suggest  a  modification  of  the  Toe- 
piitz  approximation  methoa  for  estimating  frequencies  of 
multiple  sinusoids  from  covariance  measurements.  The 
method  constructs  a  state-feedback  matrix  following  a 
low-rank  approximation  of  the  Toeoiitz  covariance  matrix 
via  singular  value  decomposition.  Ideally,  the  eigenvalues 
of  this  state-feedback  matrix  will  be  on  the  unit  circle  in 
the  complex  plane,  and  the  angles  that  they  make  with  the 
real  axis  will  be  equal  to  ihe  unknown  sinusoid  frequen¬ 
cies.  The  modification  proposed  here  exploits  this  prior 
knowledge  ot  the  modulus  of  the  eigenvalues,  and  guaran¬ 
tees  that  ev*n  :n  "he  presence  of  noise,  the  eigenvalues  of 
the  estimated  state-feedback  matrix  will  lie  on  the  unit 
circle. 


1.  INTRODUCTION 

The  problem  of  retrieving  multiple  sinusoids  (with 
frequencies  close  to  each  other)  from  perturbed  time-senes 
or  covariance  mlormation  is  of  special  mierest  in  a  vast 
range  of  signal-processing  applications.  Very  often  the 
covariance  seauence  mav  have  'o  be  estimated  from  time- 
senes  data,  as  in  Doopier  processing  in  radar.  It  is  not 
uncommon,  however,  to  encounter  applications  in  which 
the  1  time-seres)  data  are  not  measurable  while  the  covari¬ 
ance  information  is  directly  available.  Such  situations 
arise  in  astronomical  star  bearing  estimation,  inter! erence 
soectroscopy .  and  some  sensor  arrav  aoplicauons. 

In  recent  years,  there  has  been  a  great  deal  of  interest 
in  model -based  sinusoid  retrieval.  Models  convert  the 
non-iinear  problem  01  estimating  the  sinusoid  frequencies 
into  a  simoier  nroblem  01  estimating  the  parameters  of  a 
linear  moue:  ( 1 1.  The  second  steo  in  all  model-based 
methods  is  the  extraction  of  the  desired  information  (The 
frequencies i  from  'he  estimated  model  parameters  (21. 
Both  mens  are  mnor'ani  lor  ’he  overall  success  of  a 
motlel-bascc  method.  II! -conditioning  it  either  steo  can 
adversely  meet  the  overall  performance  ol  ’lie  meihod 
anu  snotud  re  lvmueu.  The  reliability  ol  the  first  steo 
Jenencis  >n  he  estimation  proceuure.  and  'hat  ol  'he 
second  steo  on  the  sensitivity  ol  the  desired  mlormation  to 
the  model  parameters  f 3 1.  A  popular  model  for  'he 
sinusoid  retrieval  nroblem  is  ihe  linear  prediction  model 
first  used  bv  Prnnv  m  i8SI. 


y(t)  -  £  a*y(t-k) 

t=i 

whose  parameters  may  be  reliably  estimated  by  the 
method  of  Tufts  ana  Kumaresan  (J|.  The  roots  of  the 
polynomial  formed  from  these  parameters  are  ideally 
expected  to  be  on  the  unit  ctrcie  in  the  comoiex  plane,  ana 
the  angles  that  they  make  with  the  real  axis  should  equal 
the  sinusoid  frequencies. 

2:  STATE-SPACE  REPRESENTATION 
It  turns  out  that  the  sinusoidal  model  is  a  very  spe¬ 
cial  case  of  the  geui.-i  linear  rational  model,  and  that  just 
as  there  are  alternate  parametenzations  of  linear  systems, 
there  also  are  alternate  parameter  sets  for  the  sinusoidal 
mode!  as  well.  Just  as  there  is  a  state-space  representation 
for  every  realization  of  a  linear,  rational  system,  there  is 
also  a  state-space  representation  for  every  realization  of 
the  sinusoidal  model.  The  state-space  representation  of  the 
special  model  for  sinusoidal  signals  (frequencies:  co;. 
i-1.2 . n)  is: 

xik  +  1)  =  Fx(k) 
y(  k )  =  hxi  k ) ' 

where  the  order  of  the  model  p  is  twice  the  number  of 
sinusoids,  and  the  eigenvalues  of  F  are  of  unit  magnitude 

and  equal  e =  "L  i- 1 .2 . a.  The  sinusoidal  signal  y(  t )  is  the 

model's  zero-input  resoonse  to  some  non-zero  initial  condi¬ 
tion  xtO).  In  I  act.  we  nave 

y(t)  =  hF'bdO).  t>0. 

and  the  covariance  rt  m)  of  the  sinusoidal  signal  satisfies 

r(m i  =  hF“"Ph:  m  AO  ( 1 ) 

where  P  is  the  state-covariance  matrix,  anu  the  superscript 
t  denotes  the  Hermitian  transoose. 

The  linear  prediction  model  is  a  canonical  realization 
of  the  above,  with 


xi  l  i  = 


F  =■ 
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ABSTRACT 

In  this  paper,  we  examine  the  performance  of  Singular 
'•  olac  Decomposition  (SVD)  based  methods  for  esumatinc  the 
frequencies  of  multiple  sinusoids.  The  concept  of  angle 
between  subspaces  is  used  to  derive  a  criteria  for  determining 
when  SV’D  based  procedures  fail.  The  signal  to  noise  ratio 
(SNR)  at  which  a  method  breaks  down  is  termed  the  threshold 
SNR  of  the  method.  It  is  then  shown  that  existing  SVD  based 
methods  have  a  higher  threshold  SNR  than  predicted  by  this 
criteria.  A  method  that  utilizes  the  singular  vectors,  and  directly 
minimizes  the  angle  between  subspaces  is  developed.  The 
method  is  shown  to  have  better  performance  at  low  signal  to 
noise  rauos.  The  procedure  lowers  the  threshold  SNR  thereby 
extending  the  range  of  SNR  for  which  SVD  can  be  used. 


INTRODUCTION 

Singular  Value  Decomposition  (SVD)  has  become  a  very 
popular  tool  in  signal  processing  problems.  In  this  paper  we 
examine  its  use  for  tne  problem  of  estimating  the  frequencies  of 
multiple  sinusoids  in  white  noise.  SVD  has  been  found  to  be  a 
very  useful  tool  in  obtaining  reliable  estimates  of  the  frequen¬ 
cies.  However,  very  little  is  known  about  certain  fundamental 
issues  relating  to  us  use.  The  issues  considered  in  this  paper  are: 
1 )  When  do  SVD  based  methods  fail? 

21  Do  the  SVD  based  methods  that  exist  for  the  sinusoid  prob¬ 
lem  make  the  best  use  of  the  decomposition? 

3)  If  not.  are  there  procedures  that  can  make  the  most  out  of 
using  SVD-7 

In  ths  paper,  a  criteria  based  on  the  angle  between  subspaces  is 
considered  to  determine  when  the  SVD  based  methods  fail. 
The  sienal  to  noise  ratio  (SNRi  ai  which  a  method  breaks  down 
is  termed  the  threshold  SNR  of  the  method.  A  low  threshold 
SNR  is  a  desirable  feature  for  a  method.  A  SVD  based  State 
Space  approach  is  considered  |9,12],  and  is  shown  to  have  a 
higher  threshold  SNR  ’ban  predicted  by  this  criteria.  The  same 
is  true  for  the  SVD  based  Linear  Prediction  method  [3. 7).  A 
method  that  utilizes  the  singular  vectors,  and  directly  minimizes 
the  angle  between  subspaces  is  then  developed.  The  method  is 
shown  to  have  better  performance  at  low  signal  to  noise  ratios, 
and  has  a  lower  threshold  SNR  The  results  also  indicate  ihai 
existing  SVD  based  methods  do  not  make  the  best  use  of  the 
subspaces. 

U.  BACKGROUND 

Here  we  briefly  describe  a  SVD  based  state  space 
approach  for  the  sinusoid  frequency  estimation  problem  19.12], 
For  this  discussion  we  assume  that  the  data  is  noise  free  and 
address  the  issue  of  noise  later.  For  the  sinusoid  problem,  the 
data  is  the  sum  of  complex  exponentials,  i.e. 

iff!  =  £  c,  e'm‘. 
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where  to,  are  the  frequencies  and  c,  are  the  complex  ampluude 
with  C,  =  Ic,  I  e'*\  ®,  being  the  phase.  Such  a  signal  can  be 

modeled  by  the  following  state  space  model 

-V, 

xl*)  =  fi  X,  , 

where  the  eigenvalues  of  F  are  of  unit  amplitude  and  — 

e'a  [12.9].  this  can  be  verified  by  considering  [he  following 
useful  realization.  6 

F  ~  riiaq  (e;“),  ()a) 

XJ>  =  <  c  i  .c  2 . 9)  and  *=(l.l I).  (lb) 

The  triple  (FX0A)  characterize  the  model  and  are  not  unique.  If 
(FA 0A)  is  one  representation  then  ( J-'FT.T-'X^hT )  is  also  a 
representation  for  any  nonsingular  matrix  r.  This 

nonuniqueness  of  the  parameter  set  enables  one  to  choose  a 
coordinate  system  that  is  less  sensitive  to  oenurbanons 
113.9,12].  The  state  space  parameters  can  be  estimated  bv  using 
the  factorization  of  matrix  D  given  below. 
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where  9  is  shown  below. 


h 

hF 
hF 2 


The  matrix  D  arises  naturally  in  a  forward  and  backward 
approach  [3].  The  matrix  F  can  be  estimated  from  9  in  the  fol¬ 
lowing  manner: 

9,  F=9j  (3) 


where 
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Perturbation  Analysis  of  an  SVD-Based  Linear 
Prediction  Method  for  Estimatina  the 
Frequencies  of  Multiple  Sinusoids 


BHASKAR  D.  RAO.  member,  ieee 


\buract — Model-based  ^oectrai  estimation  techniques  consist  of  es- 
seutiailv  two  vteps.  File  first  step  is  the  estimation  of  a  parameter  set 
from  the  data,  anti  t lie  second  step  consists  of  extracting  the  relevant 
information  trom  the  parameter  set.  A  numerical  analysis  of  (he  over¬ 
all  procedure  can  he  performed  by  conducting  a  perturbation  analysis 
o|  the  m»»  steps  separately .  We  demonstrate  this  by  study  ing  the  linear 
prediction  approach  for  estimating  the  frequencies  of  sinusoids  in  white 
noise.  It  is  shown  that  in  the  first  step,  the  continuity  of  the  generalized 
inverse  and  the  concept  of  angle  between  suhspaces  piav  an  important 
role,  (  he  continuity  concept  helps  explain  the  need  for  a  low  rank  ap¬ 
proximation.  jnd  the  quality  of  the  approximant  is  appraised  bv  using 
the  not  'on  of  angle  between  suhspaces.  For  the  second  step,  the  sensi¬ 
tive  of  the  zeros  of  ihe  predictor  polynomial  becomes  an  important 
consul*. i  atom  and  is  examined.  It  is  shown  ihal  increasing  the  order  ot 
(he  predictor  polynomial  and  computing  the  minimum  norm  solution 
pros  ides  a  mechanism  to  reduce  parameter  sensitivity. 


I.  Introduction 

IN  recent  years.  there  has  been  a  Hurry  of  activity  in  the 
applications  ot  models  to  spectral  estimation,  and  (hey 
show  great  promise  (!|.  There  is  a  multitude  of  tech¬ 
niques  with  asymptotically  similar  behavior,  and  so  their 
performance  for  tirute  data  records  becomes  an  important 
consideration.  Comparisons  based  on  simulations  have 
been  used  to  evaluate  the  methods.  In  addition  to  the  sta¬ 
tistical  behavior,  a  numerical  examination  is  useful  to  un¬ 
derstand  their  reliability/perforniaiice  in  a  finite-precision 
environment,  i.e  .  when  fixed-point  or  floating-point 
arithmetic  is  used.  In  this  paper  we  consider  this  issue. 

To  begin  'he  inaiyMs.  we  note  that  the  model-based 
methods  can.  in  general,  be  divided  into  two  distinct  steps 
as  in  Fig  1 .  Step  I  consists  of  estimating  a  parameter  set 
that  describes  the  model,  and  Step  II  consists  of  deter¬ 
minin'.:  die  relevant  information  from  the  parameter  esti¬ 
mates.  Buth  Meps  are  important  to  i he  overall  mcccss  of 
the  method.  Ill-uinuiuuning  at  either  step  can  adversely 
effect  the  overall  performance  ot  the  method  and  needs  to 
be  avoided.  Reliable  estimation  < Stop  h  is  dearly  an  im¬ 
portant  issue  and  has  been  given  a  lair  amount  of  atten¬ 
tion  For  exarr.Dle.  in  t he  deterministic  least  squares  prob- 
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Fig.  I  The  two  sieps  in  a  mouei-baseu  spectrum  estimation  method 

lem.  reliable  estimation  procedures,  e.g..  the'  QR  and 
SVD  approaches,  have  been  developed  [2],  [3|.  How¬ 
ever.  for  the  spectral  estimation  procedure,  in  addition  to 
the  tirst  step,  the  second  step  is  also  an  important  factor. 
If  the  mapping  from  the  parameter  space  to  the  spectral 
domain  is  ill-conditioned,  then  improving  Step  1  is  less 
effective.  On  the  other  hand,  a  reliable  parameter  set  may 
not  always  be  easy  to  estimate.  So  a  compromise  is  nec¬ 
essary  and  a  proper  combination  of  the  two  steps  is  im¬ 
portant  for  success.  An  evaluation  of  the  two  steps  sepa¬ 
rately  provides  a  procedure  tor  evaluating  and  comparing 
methods,  and  results  in  interesting  insignts. 

Here  we  '-(insider  the  problem  ot  estimating  frequencies 
of  sinusoids  in  white  noise.  For  the  sinusoid  problem,  let 
us  denote  the  parameter  sec  by  ff  =  [#,.  ffz.  ■  ■  •  .  0. ]  and 
the  quantity  of  interest  are  the  frequencies  w.'s.  An  error 
AO  in  the  estimate  of  0  results  in  a  error  Aw.  in  the  fre¬ 
quency  where 

r-i  l3w. 

Aw.  =  _  —  At).  —  Vw  AO 
/  =  i  off. 


i  '3w.  3w. 


Vu  =  — . 


off,  'W- 


At)  =  [Ay,,  a y>. 


The  error  At)  is  a  resuit  or  the  additive  noise  in  the  lata, 
ami  is  also  due  to  finite  precision  arithmetic 

E\  Aw. ; :  =  Cw/T!  Ay  Ay'  I  Tw  <  Am  i  /T. ,,  ' 
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ABSTRACT 

in  mis  tracer.  we  ar.ah  ze  '.he  rertormance  of  two  subsoace  based 
.memcxis.  ESPRIT.  and  the  Toeoiitz  Approximation  Method  (TAM), 
for  estimating  the  direction  ot  amvai  -DOA)  of  plane  waves  in 
•.vntte  noise,  in  the  case  of  a  linear  eqmsoaced  sensor  array,  it  is 
<nown  mat  ESPRIT,  like  TAM.  can  be  treated  via  a  state  space  for¬ 
mulation  leading  o  interesting  insights.  The  state  transition  matrix 
estimated  by  ESPRIT  is  shown  to  be  related  to  that  ootained  by 
TAM  bv  a  diagonal  similarity  transformation.  The  performance  or" 
these  mernoos  is  also  anaivzea.  In  particular,  aysmmouc  resuits  for 
the  mean  sauared  error  :n  the  estimates  of  the  direcnon  of  arrival  are 
Jenveu  for  eacn  ot  the  above  mentioned  metnotis.  Simoie  closed 
form  expressions  tor  the  one  source  case  are  derived,  and  comoared 
wttn  tne  corresponding  expressions  for  the  Minimum-Norm  metnod 
and  .MUSIC.  Comcuter  simulations  are  provided  to  substantiate  the 
analysis. 

INTRODUCTION 

Subsoace  based  metnoas  have  been  develoced  and  studied  by 
a  r.utnoer  ot  researeners.  In  this  paper  we  analvze  the  oerformance 
ot  the  sucsoace  trxsed  metr.ods  used  in  esamaang  tne  Direcnon  ot 
Amvai  iDOA)  of  plane  waves  in  noise.  Due  to  length  considera- 
nons.  in  tms  pacer,  we  present  only  a  detailed  examinanon  of 
ESPRIT  ill,  and  tne  Toeoiitz  Approximation  Metnod  (TAM)  [3]. 
The  aoproacn  can  be  extended  to  analyze  the  Minimum- Norm 
metntxi  12.10)  ana  MUSIC  1-1.5], 

MUSIC  .vas  tne  rirsr  method  that  showed  the  benefits  of  using 
a  suosnace  cased  aocroacn  141.  Some  theoreucal  resuits  companng 
MUSIC,  and  the  .Minimum-Norm  method  can  be  found  in  (6.7 j 
wnerem  a  cnaracter.zauon  of  the  methods  was  done  bv  examining 
the  null  spectrum.  Some  comparisons  of  MUSIC  with  ESPRIT  can 
be  found  :n  |(U.  ESPRIT  is.  like  MUSIC,  a  general  approach,  and 
was  deveio cea  to  overcome  some  of  the  coraoutanon  and  prior 
tntormanon  renuirements  or  MUSIC  1 1.3.9).  It  is  concemuaUy  dif¬ 
ferent  in  that  :t  caais  for  an  arrav  of  idenacai  douoiets.  i.e.  reautres 
the  arrav  -o  possess  a  disoiacement  invariance.  In  return,  it  does  not 
need  detailed  trucrmation  aoout  the  array  geometry  and  Hemem 
charactensitics.  Here  we  nmv  consider  us  penormance  in  the  con¬ 
text  of  a  linear  ecutsoacea  sensor  arrav.  The  use  of  the  Toeoiitz 
. \rrroxirnation  Metnod  ;T.\M)  for  DOA  esnmauon  was  tint  sug- 
gesteu  in  i'l.  T,\M  is  a  metnod.  like  the  .Ylinimum-Norm  method, 
mat  was  first  deveiccea  for  the  sinusoid  freouency  esamaaon  proo- 
lem  i  3.1 1 1.  A  kev  feature  of  the  method  is  that  it  is  based  on  a  state 
space  model.  Suer,  modeling  has  benents  that  were  discussed  ;n 
1  '..12.131.  ana  wii  become  clearer  as  ihe  discussion  progresses,  in 
fact  :t  will  be  mown  that  ESPRIT  can  also  be  described  under  this 
formulation  leading  to  msignts  into  us  penormance. 

The  organization  ot  the  oaoer  is  as  follows.  A  State  Space  for¬ 
mulation  is  used  to  iiescrroe  ESPRIT,  and  esraolish  the  relationsmo 
betw een  TAM  ana  ESPR.rT.  Asvmotoucal  results  tor  the  mean 
sauared  error  tn  tne  estimates  of  the  DOA  in  eacn  of  the  two  cases 
ire  derived  and  ccmcarea.  The  results  are  specialized  for  the  one 
ource  case,  and  compared  with  the  corresponding  expressions  tor 
me  Minimum- Norm  metnod  and  MUSIC.  Simutauon  results  are 
presented  ana  tr.ev  support  tne  maivsis. 
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PROBLEM  FORMULATION 

The  problems  of  estimating  the  iirecnon  of  arrival  of  \f 
incoherent  plane  waves  incident  on  a  linear  eauispaced  array  of  L 
sensors  is  considered  in  this  oaner.  For  me  ah  ooser/anon  penod 
isnapshou,  the  spaaal  samples  ot  the  signal  plus  noise  are  given  by 

Y[=\y\k)  .y?  .  ■■■  .yi"} 

'•  =  Z/7}».  I^V0*.  +N{.  (1) 

_1  =  |  I-.  1=4 

2x  d 

where  to,-  =  ~ —  stn9.-.  d  being  the  separanon  between  sensors.  \ 
the  wavelength  of  the  incident  signal,  and  9;  the  direction  of  arrival. 
The  subspace  based  meutoas  esumate  the  signal  zeros 

*  =e;U|  ,  i  =  l . Vf.  from  whicn  tne  signal  freaucnctes  w.-'i  and  then 

the  DOA's  9;  are  determined.  As  in  (6,7],  it  is  assumed  that  is  a 
seauence  of  independent,  mean  zero  Gaussian  random  vectors,  i.e 
£[,Vt  iVf  ]  =  a;  /  by.  The  noise  is  assumed  uideoendent  of  ihe  com¬ 
plex  signal  amplitudes  of1  wnic.i  are  also  modeled  as  oeing  jointly 
Gaussian.  The  covariance  matrix  P  of  the  amplitudes  whose  ele¬ 
ments  are  P,,,  where  P;,  =  Elpf' p]'*%  is  assumed  to  be  of  rank  M 
and  has  distinct  eigenvalues. 

Now  we  make  some  comments  regarding  the  notation  and  con¬ 
ventions  used  Throughout  the  the  caper.  '7*  is  used  to  denote  tran¬ 
spose.  "  to  denote  complex  comugate.  H'  :o  denote  complex  conju¬ 
gate  transpose  and  +  to  denote  me  Moor  e-Penrose  pseudo-inverse. 
Also  £(.]  and  the  superscript  will  be  'used  mtercnangeablv  to 
denote  the  expectation  operator.  )n  this  paper,  is  used  to  denote 
estimates,  and  suosenot  s  to  denote  parameters  associated  with  the 
signal  alone.  Also  if  y  is  a  vector  of  lengtn  L.  then  y'  and  v"  denote 
vectors  of  length  L-l  formed  from  the  rirsi  and  last  L~i  elements 

of  y,  i.e.  if  i  •  >'; . >'£.-!  •  'Sl  ■'  •  -  l*crl 

y  =  [  V[  ,  y; . PL-l  i  •  Old  =  (  >'J . .  >'£  J  ■  F°r 

compactness  of  representauon.  we  derine  two  i L-ifcL  matrices  W 
ana  \V" 

(C  =  jf'c-iuiL-n  •  ®  j  W  ~  9  • /<p-ni(£.-ii  |  ■  1-^ 

Note  that  v'  =  W  y  ana  v"  =  IF"  For  matnees  we  have  the  follow¬ 
ing  convention.  If  3  is  a  LxM  matrix,  then  3  (  and  3  ->  are  <L  -!)xM 
matrices  containing  the  tint  and  last  L-l  rows  of  3.  i.e. 

r*.  l 


by 

3g 

Also  Ax  is  used  to  denote  the  error  in  uuanutv  r.  where  s  mav  be  a 
scalar,  vector  or  matnx. 

The  covanance  matnx  of  me  observation  vector  lias  an 
eigenaecomoosition  and  can  be  wntten  is 

R  =  Ft  i'V  =  Z  SiS"  ~  E  A  =  £'  *■  a':-  ' 41 
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ABSTRACT 

In  this  paper,  we  analyze  the  performance  of  Root-Music,  a 
variation  of  the  popular  MUSIC  algorithm,  for  estimating  the 
direction  of  arrival  (DOA)  of  plane  waves  in  white  noise  in 
the  case  of  a  linear  equispaced  sensor  array.  The  performance 
of  the  method  is  analyzed  by  examining  the  perturbation  in 
the  roots  of  the  polynomial  formed  in  the  intermediate  step  of 
Root-Music.  In  particular,  avsmptotic  results  for  the  mean 
squared  error  in  the  estimates  of  the  direction  of  arrival  are 
derived.  Simple  closed  form  expressions  are  derived  for  the 
one  and  two  source  case  to  get  further  insight.  Computer 
simulations  are  provided  to  substantiate  the  analysis.  An  im¬ 
portant  outcome  of  this  analysis  is  an  explanation  as  to  why 
Root-Music  is  superior  to  the  popular  MUSIC  algorithm 
where  one  examine  the  peaks  of  the  spatial  spectrum. 

L  INTRODUCTION 

Eigen-Decomposition  based  methods  have  recently  been 
extensively  used  in  estimating  the  Direction  of  Arrival  (DOA) 
of  plane  waves  in  noise  (1.6).  These  methods,  often  referred 
to  as  subspace  based  methods,  have  been  shown  to  perform 
very  well  and  are  capable  of  resolving  closely  spaced  sources. 
MUSIC  was  the  first  method  that  showed  the  benefits  of  using 
a  subspace  based  approach  [Ij.  The  MUSIC  algorithm  com¬ 
putes  a  spatial  spectrum  from  the  noise  subspace,  and  deter¬ 
mines  the  DOA's  from  the  dominant  peaks  in  the  spectrum. 
Another  popular  variation  of  MUSIC  is  Root-Music  ( 2). 
Root-Music,  as  described  in  more  detail  later,  is  similar  to 
MUSIC  in  many  respects  except  that  the  DOA's  are  deter¬ 
mined  from  the  roots  of  a  polynomial  formed  from  the  noise 
subspace.  Though  MUSIC  is  applicable  to  general  known 
array  configurations  Root-Music  is  mainly  suitable  in  the 
context  of  a  linear  equispaced  sensor  array  (2).  Some  theoret¬ 
ical  results  comparing  MUSIC,  and  the  Minimum-Norm 
method  can  be  found  in  (3.4)  wherein  a  characterization  of 
the  methods  was  done  by  examining  the  null  spectrum.  Our 
work  examines  Root-Music,  and  characterizes  the  mean 
squared  error  in  the  estimates  of  the  DOA's  directly.  The 
analysis  provides  insight  into  why  Root-Music  is  superior  to 
the  popular  MUSIC  algorithm  where  one  examines  the  peaks 
in  the  spatial  spectrum. 


II.  PROBLEM  FORMULATION 
The  problem  of  esi!iiiz',ng  the  direction  of  arrival  of  M 
incoherent  plane  waves  incident  on  a  linear  equispaced  array 
of  I.  sensors  is  considered  in  this  paper  For  the  kih  observa- 
non  period  I  snapshot),  the  spatial  samples  of  the  signal  plus 
noise  are  given  by 

V  p/*’  e-"°'  .  Z  p/*’  CJ,L'])W  *N[.  (I) 
=  1  .  =! 

’hi  *  wcrr  r>  'uopnrf,9'j  by  ARMY  ^.'^arch  Qffic? 
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where  0);  =  — - —  sinQ, ,  d  being  the  separation  between  sen* 

sors.  X  the  wavelength  of  the  incident  signal,  and  9,  the  direc¬ 
tion  of  arrival.  Root-Music  estimates  e  .  the 

signal  zeros,  from  which  w,'s.  the  signal  frequencies,  and 
then  the  DOA's  9;  are  determined.  As  in  (3.41,  the  noise 
vector  JVt  is  assumed  to  be  a  zero  mean,  complex  while 
Gaussian  random  vector,  i.e.  Nt  Nj  =  a*  /  5U.  The  noise  is 
assumed  independent  of  the  complex  signal  amplitudes  p/*1 
which  are  also  modeled  as  being  jointly  Gaussian.  The 
covariance  matrix  P  of  the  amplitudes  whose  elements  are 
P,,,  where  Pi;  =  E  ( p/4  'p," lt  ‘l.  is  assumed  to  be  of  rank  M 
and  has  distinct  eigenvalues. 

In  this  paper.  'T'  is  used  to  denote  transpose,  *  to  denote 
complex  conjugate,  and  H  to  denote  complex  conjugate  tran¬ 
spose.  Also  £[.]  and  the  superscript  "  will  be  used  inter¬ 
changeably  to  denote  the  expectation  operator.  In  this  paper. 

is  used  to  denote  estimates,  and  subscript  s  and  n  to  denoie 
parameters  associated  with  the  signal  and  noise  respectively. 

The  covariance  matrix  of  the  observation  vector  Yt  has 
an  ciccndccomposiuon  and  can  be  written  as 


R  =  Yt  Y>J  =  V  X, 

/  =  1 

S,  5/'  =  £  A  E" 

=  £,  A,  £/'  *  o-l. 

1 2) 

where 

£=(£,.  5;  . 

•  5J.  £,=(£, 

.5;  ....  5W] 

i  3a) 

A  =  dtag  (X| ,  X- ,  • 

■  .  XL ).  and 

Aj  =  ding  (Xf .  Xf . 

•  -Ky 

■  ?h) 

Also 

*•!  =  Aj'  +  c„:  >  X,  =  X;*  +  o*  >  >  ltf  =  X,’f  *  a; 

>  Vvm  =  ■  =  al  =  o;. 

X,  are  the  eigenvalues  of  R  .  and  5,  the  corresponding  ortho- 
normai  eigenvectors.  The  subspace  based  methods  esnmaie 
the  DOA’s  from  either  the  subspace  E,,  t.e.  the  space  spanned 
by  the  eigenvectors  corresponding  to  the  \1  dominant  ei sen- 
values.  termed  the  signal  subspace,  or  ns  orthogonal  comple¬ 
ment  £v  where  £v  =  (£i^|  ,  ■  ,5;  !,  termed  the  noise 

subspace.  Pj  and  Pv  arc  used  to  denote  the  orthogonal  pro¬ 
jection  operators  onto  the  signal  and  noise  subspace  respec¬ 
tively. 

This  paper  considers  the  effect  of  using  an  estimated 
covariance  matrix.  Usually  an  esumate  of  the  covariance 
matrix  is  obtained  by  (time)  averaging  N  independent 
snapshots,  i.e. 

R  =  —  V  Yt  Y ?  =  £  A  E11, 

■v  *7i 

where 

£  =  (5,  .  5j  .  •  ■  •  .  SL  ],  A  =  ding  (Xt ) 
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ABSTRACT 

In  Bus  paper,  we  examine  state  space  methods  for  estimating  the  direction  of  arrival  of  plane  waves  from  a 
linear  equispaced  sensor  array.  From  a  computational  point  of  view,  the  approach  only  involves  matrix  operations  and 
is  suitable  for  systoiicwavetronc  impiementauon. 

I.  INTRO DL'CTION 

Development  of  signal  processing  algorithms  for  implementation  on  Systolic/'vVavetrcnt  architecrunes  have  been 
of  interest  in  the  past  few  years  [IRC].  In  this  pacer,  we  examine  methods  for  estimating  the  direction  of  arrival 
(DOA)  of  plane  waves  using  a  Linear  equispaced  sensor  array  with  a  view  to  implement  them  on  svstoiic/wavetront 
arrays.  In  the  context  of  DOA  esnmaben.  subspace  based  methods  have  recently  received  a  great  deal  of  attention. 
Examples  of  such  methods  are  MUSIC  [41,  ESPRIT  [51,  the  Taeplitz  Approximauon  Method  (TAM)  [6]  ex.  In  par- 
Gcular  MUSIC  has  received  a  great  deai  of  attention.  The  various  steps  of  a  MUSIC  algorithm  are  summarized  in 
Fig.  I.  Note  that  there  are  efficient  systolic  architectures  to  perform  the  first  two  steps,  namely  covariance  computa¬ 
tion  and  eigendecomposiucn.  However  the  process  of  computing  the  null  spectrum  D  (to)  and  locating  the  peaks  in 
its  inverse  is  computationally  expensive  and  not  very  amenabie  to  systolic  processing.  In  this  paper,  we  examine 
state  space,  based  methods.  State  space  methods  for  these  problems  were  first  suggested  in  [6.7].  The  computauonal 
steps  m  the  method  are  based  soieiy  on  matrix  operauons  making  it  an  attractive  candidate  for  systolic  implementa¬ 
tions.  Note  that  ESPRIT  for  a  linear  equispaced  array  can  be  treated  via  this  formalism.  Because  of  its  applicability 
to  more  general  army  configurauons,  the  ideas  presented  here  have  'wider  applicability. 


O.  STATE  SPACE  APPROACH 

The  problem  of  estimating  the  direction  of  arrival  of  \t  incoherent  plane  waves  incident  on  a  linear  equispaced 
array  of  L  sensors  using  sxte  space  methods  is  considered  in  this  paper.  For  the  <:h  observation  period  (snapsnot), 
the  spatial  samples  of  the  signal  plus  noise  are  given  by 
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abstract 

In  rhis  paper,  the  roundoff  noise  properties  oi  floating 
point  statc-sDace  digitai  niters  is  examined.  It  is  shown  that 
the  roundoff  notse  behavior  of  these  niters  is  related  to  their 
coefficient  sensitivity.  Expression  for  the  variance  of  the  out¬ 
put  roundoff  noise  is  aenved  when  the  input  to  the  niter  is  a 
zero  mean  wide  sense  stationary  random  process.  These 
expressions,  as  expected,  depend  on  the  input  signal  statis¬ 
tics  along  with  tne  niters  parameters,  in  particular  the  con- 
troiiaoiiity  and  ooservabiiitv  grammians. 


I.  INTRODUCTION 

Floating  point  arithmetic  is  widely  used  in  signal  pro¬ 
cessing  computations.  More  recently,  with  the  increasing 
availability  of  floating  point  capability  m  signal  processing 
cnips.  msignt  into  algorithms  employing  floating  point  arith¬ 
metic  is  of  interest.  This  paper  considers  the  probiem  of 
digital  filter  implementation,  and  examines  the  effect  of 
floating  point  arithmetic  on  them,  [n  particular  floaang  point 
state-space  digitai  filters  are  analyzed.  Tne  prooiem  of  using 
axed  point  arithmetic  and  design  of  minimal  roundoff  filters 
;s  weii  known  [1,21.  in  contrast,  little  is  mown  aoout 
rounc  off  noise  in  filters  employing  fioanng  point  antnmeuc. 
Effect  of  floating  point  antnmenc  on  digital  niters  has  been 
studied  in  tne  context  of  direct,  cascade  and  parallel  forms 
[3,-i],  and  parnaily  for  state  space  digital  niters  [5].  Here  a 
detailed  treatment  of  tne  effects  of  using  floating  point  arith¬ 
metic  in  state-space  digitai  filters  is  presented. 

II.  floating  point  arithmetic 

Throughout  this  paper  it  will  be  assumed  that  floating 
point  numoers  are  stored  in  the  form  (sign  mi. 2'.  where  a 
and  v  have  a  fixed  numoer  of  bits.  Also  it  wiil  be  assumed 
that  rounding  is  used  in  all  operations.  Truncation  can  be 
deal;  with  in  a  similar  manner.  Tne  notation  /if.,  will  be 
used  :o  denote  the  machine  number  resulting  from  fioanng 
point  operation,  in  fioanng  point  operations  There  are  errors 
in  noth  additions  and  multiplications  [3,4],  i.e. 


and 

/.Ary  '«!)•'!  -01. 

a^d  -2”  Se;.5<2'v,  g  being  tne  number  of  bits  used  to 
represent  tne  mantissa,  the  e.ror  vanaoies  e  and  0  are 
assumed  to  oe  rancom  variables  uniformly  distributed 
oetween  -2"*  and  IT  The  error  vanaoies  have  zero  mean 

‘This  worx  was  supported  by  -he  ARMY  Researc.n  Office 
under  Grant  No.  DAAL-03-86-K.-0107.  and  by  the  National 
Science  Foundation  under  Grant  No.  MIP-S7- 1 I9Sd. 


and  vanance  a."  =  — .  Tr.e  error  :n  the  case  of  addition  is 

bx  ->■)  and  m  the  case  of  multiplication  is  or  v,  A 
significant  difference  between  fixed  point  and  floating  point 
arithmetic  is  tha'  me  error  caused  by  rouncina  in  floating 
point  arithmetic  is  dependent  on  the  signal  Tins  depen¬ 
dence  of  the  rounding  error  on  the  signal  makes  the  anaivsis 
less  Uactaoie.  To  cany  out  the  analysts  some  carefui  mam- 
pulanons  are  necessary.  We  now  consider  in  detail  the  float¬ 
ing  point  computation  of  the  inner  product,  a  computation 
often  used  in  a  state-space  filter.  The  inner  product  of  tne 
vector  j  =  ,u.  .  j-..i5]'  and  x  =  [xt  .jtj.xjJ'  is  given  bv 

>•  =  ±:  *  *  3,X.  - J:X;  -  JjX,. 

As  will  be  clear  irom  the  discussion  'oeiow.  the  error  due  to 
n.nite  precision  floating  point  arithmetic  will  depend  cn  tne 
order  of  the  operations. 

y :  =  /i(/ifa ,x ~;7 (/ifu-x-,  -  /!(u,xy> ;  1  Gai 

=  '0 -x  j(l  —  5.,  -  u -r .. :  -  i-,  -  j ,x  ji 1  -  -Si j!(1  -  e.  1)  ,  i  -  i-j 

u ; i »  —  A, ,  a; .  -*-  J  -t  I  —  A- ,  J  ^  —  Jiil  —  A-^ , x i . 

=  ,j.  -  A Uiiix,  -  :a  -  -  Aj^X;  -  ij,  -  Ad  1.1X3 

=  ■  d  -  Aa  /  X  .  :  t  b  1 

where  Au:  =  A,  j,  and  A,  =  5;  -  and  A;  =  5-  -  s,  -  £;  and 
Ai  =  -  e;  Sucn  a  first  order  approximation  ;s  reason¬ 

ably  accurate  since  the  errors  will  be  assumed  to  be  d. 
random  variables.  Such  approximations  will  be  mace 
throuen  this  paper.  Note  that*  ,.Au  ,*  =  2j:pT. 
(Ac;,-  =  3 j-  or  _  and  Auy~  =  3  j,:  or.  Also 

Au .  aj  -  =  j  ■  j  -  or.  Au  ■  Au ,  -  j  ■  j ,  err.  and 

Au:AUi  =  2 Note  dial  tne  resulting  errors  are  corre¬ 
lated.  If  the  operations  of  the  inner  procuc:  were  done  m  a 
different  order,  men  tne  resu.tant  errors  have  different 
cnaractensttcs.  i.e. 

fl<  fi'a-x-i  -  hij-Oi  -  *  /;,  jj,, 

=  a-  *  a'j,,x;  -  ;j.  •  A'j.ii,  -  a,  -  a'jmt. 

1  he  errors  Aj  ire  correlated  and  have  different  statistics. 
An  .mrenant  ooservanen  in  this  context  is  mat  tne  inner 
procuc:  obtained  using  floating  point  amtnmetic  is  tne  -exact 
inner  product  of  a  perraroed  vector  i.e.  j  -  aj.  me  In 
this  paper  we  wiii  snow-  that  consioeraoie  tnsignt  can  re 
obtained  py  using  me  aoove  simpie  interpretation  of  me 
floating  point  inner  product  computation  process.  Note  mat 
the  statistics  ot  aj  depend  on  the  order  n  -.vnicn  me  m.r.er 
product  is  comDutec.  W'e  now  appiy  this  mode;  to  me  proo¬ 
iem  of  state  space  digital  niters. 


Tne  overoar  denotes  the  expectation  operator 
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ABSTRACT 

In  this  paper,  we  analyze  the  performance  of  the 
Minimum-Norm  method  for  estimating  the  airecaon  of  arrival 
CDOA)  of  piane  waves  in  white  noise  in  the  case  of  a  iinear 
cquispaced  sensor  array.  The  performance  of  the  method  is 
anaivzed  by  examining  the  perturbation  in  the  roots  of  the 
polynomial  formed  in  the  intermediate  step  of  the  Minimum- 
Norm  method.  In  particular,  avsmptonc  resuits  for  the  mean 
squared  error  in  the  estimates  of  the  direction  of  arrival  are 
derived.  Simple  closed  form  expressions  are  derived  for  the 
one  and  two  source  case  to  get  further  msignt.  Computer 
simulations  are  provided  to  substantiate  the  analysis.  The 
resuits  obtained  are  compared  to  those  obtained  for  Root- 
Music  and  it  is  shown  that  the  relative  performance  of  the  two 
methods  is  directly  dependent  on  tne  ratio  of  their  parameter 
sensitivities. 

L  INTRODUCTION 

Eigen-Decomposition  based  methods  have  recently  been 
extensively  used  in  estimating  the  Direction  of  .Arrival  fDOA) 
of  piane  waves  in  noise.  These  methods,  often  referred  to  as 
subspace  based  methods,  have  been  shown  to  perform  very 
well,  and  are  capable  of  resolving  closely  spaced  sources.  In 
recent  years  a  statistical  evaluation  of  these  methods  has  been 
conducted  by  a  number  of  researchers  [2-4.9].  Our  work 
examines  the  Minimum-Norm  method  [I],  and  characterizes 
the  mean  squared  error  in  the  estimates  of  the  DOA's. 
Motivated  by  our  observations  regarding  Root-MUSIC  [4],  we 
examine  the  error  in  the  roots  of  the  polynomial  formed  in  the 
intermediate  step  of  the  Minimum-Norm  method.  The  results 
obtained  are  compared  with  those  for  Root-Music  leading  to 
interesting  insights. 

0.  PROBLEM  FORMULATION 

The  problem  of  esumaang  the  direction  of  arrival  of  M 
incoherent  plane  waves  incident  on  a  iinear  tauispaced  array 
of  L  sensors  is  considered  in  this  paper.  For  the  <tn  observa¬ 
tion  period  (snapshot),  the  spatial  samples  of  the  signal  plus 
noise  are  given  by 

r  v  m  w  ..  ,,  i 

Y\  -  \  ZPi‘ )-Zp!‘ -,Ykr.  (1) 

h  =  l  i=l  (  =  1 

where  CO,  =  sm0; .  d  being  the  separation  between  sen- 

A. 

sors.  X  the  waveiength  of  the  incident  signal,  and  0,  the  direc¬ 
tion  of  arrival.  As  in  [2.3],  the  noise  vector  is  assumed  to 
be  i  zero  mean,  complex  white  Gaussian  random  vector,  i.e. 
•V<  N;  =  a‘  /  ou.  The  noise  is  assumed  independent  of  the 
complex  signal  amplitudes  p,'k)  which  are  also  modeled  as 
being  jointly  Gaussian.  The  covariance  matrix  °  of  the 
amplitudes  whose  elements  are  P„,  where  P:/  = 
is  assumed  to  be  of  rank  M  and  has  distinct  eigenvalues.  In 
this  paper,  the  overbar  wiil  be  used  to  denote  die  expec- 
tauon  operator* 

Subspace  based  methods  estimate  r,  =  e;n  .  i  =  !,..„lf , 
the  signal  zeros,  from  which  W.  ’s .  the  signal  frequencies,  and 
then  the  DOA's  0,  are  determined.  They  utilize  the  eigen- 
decomposition  of  the  covariance  matrix  of  the  observation 
vector  Yk,  i.e. 

*  T-m  s  -orK  suororteh  by  tbe  ARMY  Research  3**109 
jnce-  Grant  No.  3AAL-  33-  35-  <  -  31  37  . 


wnere  co,  = 


=  yTW  =  I  x,  S,  Si1  =  £  A  Eh  =  E,  Aj  E“  -r  z;l .  (2) 


'  •Aif). 


E  -  [5j  ,  5;  .  •  ,  $l  ].  Et  =  [S  i  .  S; ....  S\,t  ]  (3a i 

A  =  diag  (Xj ,  X^ ,  ,  Xr_).  and 

A  i  =  diag  0.  f,X£,  ■  ■  ■  ,X£).  (3b) 

Also 

M  =  Xf  -  a~  >  -  X!  ~  a*  >  ..  >  Xw  =  ~  a; 

>  =  •  •  •  =  XL  =  o-. 

X,  are  the  eigenvalues  of  R ,  and  S:  the  corresponding  ortho- 
normal  eigenvectors.  This  paper  considers  the  effect  of  using 
an  estimated  covanance  matrix.  Usually  an  esumate  of  the 
covariance  matrix  ts  obtained  by  (tame)  avenging  :V  indepen¬ 
dent  snapshots,  Le. 

1  <v  .  . 

R  =  —  V  r  Y?  =  £  A  £  ", 

N  i=i 
where 

£  =  [  S  j  .  Sj  ,  '  ’  ‘  ,  S[_  ],  A  =  diag  (Xt ). 

Let  St  -  St  +  and  £.t  =  Xk  +  The  analysis  makes  use 
of  the  asymptotic  properties  of  the  errors  i"it  and  3*  ■  Asymp¬ 
totically  the  errors  are  jointly  Gaussian  nndom  variables  with 
the  errors  in  the  eigenvectors  being  independent  of  the  errors 
in  the  eigenvalues  [5].  It  has  been  shown  in  [5]  that  for  the 
errors  corresponding  to  the  M  signal  eigenvectors, 

TU  X  ■  -  Sr  S'J  5U  +  a (,V-1 ).  (4) 

N  r  =  l  (Xt  -  A.,)‘ 


Tun/=-  — — — t~t  S/  sk  (i  -  5^ )  -r-  o  ( ,v~l ),  (5 

N  (Xt  -  X,)* 

where  is  the  kronecker  delta.  In  [2.3],  it  was  shown  that 
— -  A  X, 


1V  t~i  (X*  -x,r 


St  =  a*  Sk  *  o 


HI.  MINIMUM-NORM  METHOD 
In  the  Minimum-Norm  method  [1],  the  DOA's  are  found 
by  iocanng  the  peaks  of  S( co),  where 

S(w)  -  — — — , 

D  (co) 

and  D( col.  termed  the  null  spectrum  [2  3],  is  defined  as 
D  i cm  =  ![l.gw]  V'tcoll2. 

with  VCcol  =  —=r  (1  ,  eja>  ,  •  •  •  ,  ,>/ (4 -'•>“]*  [1  ,  g*j*  is  a 

~  -  aenotes  the  pseudoinverse,  T  denotes  transpose.  * 

denotes  compiex  conjugate.  H  denotes  compiex  conjueate 
transpose.  Aiso  *  is  used  to  denote  estimates  and  subscript  s 
and  n  denote  parameters  associated  with  signal  and  noise 

resoectiveiv. 
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SENSITIVITY  CONSIDERATIONS  IN  STATE  SPACE  MODEL  t 
BASED  HARMONIC  RETRIEVAL  METHODS 


Bhaskar  D.  Rao 
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ABSTRACT 

In  tiiis  paper,  state  space  models  are  considered  for  estimating  the  frequencies  of  multipie  sin- 
susoids,  and  robust  coordinate  systems  for  frequency  estimation  are  identified.  It  is  shown  that 
the  ideal  parameter  set  involves  estimadng  an  unitary'  state  transition  matrix,  and  then  comput¬ 
ing  its  eigenvalues  to  obtain  an  estimate  of  the  frequencies.  Procedures  to  estimate  such  ma¬ 
trices  from  covariance  and  time  series  data  are  examined.  It  is  shown  that  two  state  space 
methods,  the  Toeplitz  Approximation  Method  (TAM)  and  the  Direct  Data  Approximation 
(DDA),  estimate  robust  state  transition  matrices. 


This  work  was  supported  bv  the  ARMY  Research  Office  under  Grant  No. 
DAAL-03-86-K-0107. 
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PERFORMANCE  ANALYSIS  OF  ESPRIT  AND  TAM  IN  DETERMINING 
THE  DIRECTION  OF  ARRIVAL  OF  PLANE  WAVES  IN  NOISE 
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ABSTRACT 

In  this  paper,  two  subspace  based  methods.  ESPRIT  and  the  Toeplitz  Approximation 
Method  (TAM),  for  estimating  the  direction  of  arrival  (DOA)  of  plane  waves  in  white  noise 
in  the  case  of  a  linear  equispaced  sensor  array  are  evaluated.  It  is  shown  that  the  least  squares 
version  of  ESPRIT  and  TAM  result  in  the  same  estimate,  and  are  statistically  equivalent.  It  is 
shown  that  asymptotically,  the  estimates  obtained  using  Least  Squares  ESPRIT  and  Total 
Least  Squares  ESPRIT  have  the  same  mean  squared  error.  Expressions  for  the  asymptotic 
mean  squared  error  in  the  estimates  of  the  direction  of  arrival  are  derived  for  the  methods. 
Simple  closed  form  expressions  are  derived  for  the  one  and  two  source  case  to  get  further 
insight.  Computer  simulations  are  provided  to  substantiate  the  analysis. 


Hus  work  was  supported  hv  the  ARMY  Research  Office  under  Grant  No.  DAAL-03-86-K- 
0107. 
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ABSTRACT 

In  this  paper,  we  analyze  the  performance  of  Root-Music,  a  variation  of  the  popular 
MUSIC  algorithm,  for  estimating  the  direction  of  arrival  (DO A)  of  plane  waves  in 
white  noise  in  the  case  of  a  linear  equispaced  sensor  array.  The  performance  of  the 
method  is  analyzed  by  examining  the  perturbation  in  the  roots  of  the  polynomial 
formed  in  the  intermediate  step  of  Root-Music.  In  particular,  avsmptouc  results  for  the 
mean  squared  error  in  the  estimates  of  the  direction  of  arrival  are  derived.  Simpie 
closed  form  expressions  are  derived  for  the  one  and  two  source  case  to  get  further  in¬ 
sight.  Computer  simulations  are  provided  to  substantiate  the  analysis.  An  important 
outcome  of  this  analysis  is  an  explanation  as  to  why  Root-Music  is  superior  to  the  po¬ 
pular  MUSIC  algorithm  where  one  examines  the  peaks  of  the  spatial  spectrum. 
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ROUNDOFF  NOISE  IN  FLOATING  POINT  STATE  SPACE  DIGITAL  FILTERS* 
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ABSTRACT 

In  this  paper,  the  roundoff  noise  properties  of  floating  point  state-space  digital  filters  is 

examined.  It  is  shown  that  the  roundoff  noise  behavior  of  these  filters  is  related  to  their 

coefficient  sensitivity.  An  exact  expression  for  the  variance  of  the  output  roundoff  noise  is 
-  » 

derived  when  the  input  to  the  filter  is  a  zero  mean  wide  sense  stationary  random  process. 
Surprisingly,  these  expressions  are  tractable  and  depend  on  the  input  signal  statistics  along 
with  the  filter  parameters,  in  particular  the  controllability  and  observability  grammians.  For  the 
case  where  double  precision  accumulation  is  used,  it  is  shown  that  the  optimal  filters  are  simi¬ 
lar  to  those  obtained  when  using  fixed  point  arithmetic.  Also  the  roundoff  noise  of  (single  pre¬ 
cision)  second  order  filter  structures  is  studied  in  detail. 
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Abstract:  The  performance  of  the  minimum-norm 
method,  both  root  and  spectral  forms,  for  estimat¬ 
ing  the  direction  of  arrival  (DOA)  of  plane  waves 
in  white  noise  is  analysed  for  the  case  of  a  linear 
equispaced  sensor  array.  In  particular,  asymptotic 
results  for  the  mean  squared  error  in  the  estimates 
of  the  direction  of  arrival  are  derived.  Simple 
closed-form  expressions  are  derived  for  the  one 
and  two  source  case  to  get  further  insight.  Com¬ 
puter  simulations  are  provided  to  substantiate  the 
analysis.  The  results  ootained  are  compared  to 
those  obtained  for  root-MUSIC  and  it  is  shown 
that  the  relative  performance  of  the  two  methods 
is  directly  dependent  on  the  ratio  of  their  param¬ 
eter  sensitivities. 


1  Introduction 

In  this  paper,  we  analyse  the  statistical  performance  of  an 
eigendecomposition  based  method,  the  minimum-norm 
method  [1],  Eigendecomposition-based  methods  have 
recently  been  extensively  used  in  estimating  the  direction 
of  arrival  (DOA)  of  plane  waves  in  noise.  These  methods, 
often  referred  to  as  subspace-based  methods,  have  been 
shown  to  perform  very  well,  and  are  capable  of  resolving 
closely  spaced  sources.  In  recent  years  a  statistical  evalu¬ 
ation  of  these  methods  has  been  conducted  by  a  number 
of  researchers  [2-6],  For  instance,  some  theoretical 
results  comparing  MUSIC,  and  the  minimum-norm 
method  can  be  found  in  References  2  and  3,  wherein  a 
characterisation  of  the  methods  was  done  by  examining 
the  null  spectrum.  More  recent  work  on  the  analysis  of 
MUSIC  can  be  found  in  References  4,  5  and  6.  Some 
comparisons  of  MUSIC  with  ESPRIT,  as  well  as  root- 
MUSIC  with  ESPRIT,  based  on  computer  simulations, 
can  be  found  in  References  7  and  8.  Our  work  examines 
the  minimum-norm  method,  and  characterises  the  mean 
squared  error  in  the  estimates  of  the  DOAs.  Motivated 
by  our  observations  regarding  MUSIC  [6],  we  examine 
the  error  in  the  roots  of  the  polynomial  formed  in  the 
intermediate  step  of  the  minimum-norm  method.  The 
results  obtained  are  compared  with  those  for  root- 
MUSIC,  leading  to  interesting  insights.  It  is  shown  that 
the  relative  performance  of  the  two  methods  is  directly 
dependent  on  the  ratio  of  the  parameter  sensitivities.  It  is 
shown  that  this  error  is  also  related  to  the  original 
minimum-norm  method,  in  which  one  determines  the 
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direction  of  arrival  from  the  peaks  of  a  spectrum  [1].  In 
addition,  the  analysis  gives  an  insight  into  the  difference 
between  procedures  that  employ  a  root-finding  approach 
and  those  that  examine  peaks  in  the  spectrum. 

2  Problem  formulation 

The  problem  of  estimating  the  direction  of  arrival  of  .V/ 
incoherent  plane  waves  incident  on  a  linear  equispaced 
array  of  L  sensors  is  considered  in  this  paper.  For  the  feth 
observation  period  (snapshot),  the  spatial  samples  of  the 
signal  plus  noise  are  given  by 


^  =  [yf.yf,  ...,y£’] 

-  ^ 

_i“l  i=  1 


V  ‘ 


+  SJ  (1) 


where  , 

•2 ltd  . 

cu;  =  — —  sin  0; 

A 

i i  being  the  separation  between  sensors,  A  the  wavelength 
of  the  incident  signal,  and  0,  the  direction  of  arrival. 
Subspace-based  methods  estimate  the  signal  zeros  = 
ejm,  i  =  1,  ....  M,  from  which  the  signal  frequencies  cu,- 
and  then  the  DOAs  0f  are  determined.  As  in  References  2 
to  6,  the  noise  vector  St  is  assumed  to  be  a  zero  mean, 
complex  white  Gaussian  random  vector,  i.e.  iVj.V}1  = 
a;  ISkl.  The  noise  is  assumed  to  be  independent  of  the 
complex  signal  amplitudes  p|kl  which  are  also  modelled  as 
being  jointly  Gaussian.  The  covariance  matrix  P  of  the 
amplitudes  whose  elements  are  Pijt  where  PiS  = 
[pSklp*,k>],  is  assumed  to  be  of  rank  .Vf  and  has  distinct 
eigenvalues. 

The  overbar  * — ’  will  be  used  to  denote  the  expectation 
operator  in  this  paper.  Also,  T  is  used  to  denote  trans¬ 
pose,  *  to  denote  complex  conjugate.  H  to  denote 
complex  conjugate  transpose,  and  +  to  denote  the 
Moore-Penrose  pseudo-inverse.  Further,  '  is  used  to 
denote  estimates  and  subscripts  s  and  n  to  denote  param¬ 
eters  associated  with  the  signal  and  noise  respectively.  \x 
is  used  to  denote  the  error  in  parameter  .x,  which  may  be 
a  scalar,  vector  or  matrix. 

The  covariance  matrix  of  the  observation  vector  Yk  is 
given  by 

R  =  yk  T*H  =  Vs  PV*  +•  a;  I  (2a) 

where 
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ABSTRACT 

In  this  paper,  we  study  the  relationship  between  a  state  space  approach  and  the  matrix 
pencil  method  for  the  problem  of  estimating  the  parameters  of  damped  and  undamped 
sinusoids.  It  is  shown  that  the  methods  are  very  similar  in  the  exact  data  case,  and  have  minor 
difference  in  implementation  in  the  noisy  data  case. 
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ABSTRACT 

In  this  paper,  we  analyze  the  effect  of  using  a  spatially  smoothed  forward-backward 
covariance  matrix  on  the  performance  of  Root-Music  and  the  root  Minimum-Norm  method  for 
estimating  the  direction  of  arrival  of  plane  waves  in  white  noise  in  the  case  of  a  linear 
equispaced  sensor  array.  In  particular,  avsmptotic  results  for  the  mean  squared  error  in  the 
estimates  of  the  signal  zeros  and  the  direction  of  arrival  are  derived.  Simple  approximations 
are  made  to  obtain  insight  into  the  performance  of  the  methods.  It  is  shown  that  the  use  of  a 
forward- backward  covariance  estimate  and  spatial  smoothing  is  far  more  beneficial  to  the 
Minimum-Norm  method  than  to  MUSIC.  Proper  spatial  smoothing  enables  the  performance  of 
the  Minimum-Norm  method  to  be  made  comparable  to  MUSIC.  Also,  in  the  case  of  MUSIC 
tnough  spatial  smoothing  improves  the  spectral  efficiency  factor,  minimal  smoothing  is  desir¬ 
able  for  Root-Music.  Computer  simulations  are  provided  to  substantiate  the  analysis. 
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ABSTRACT 

This  paper  reviews  state  space  model-based  methods  for  signal  processing  applications.  A  state  space  frame¬ 
work  is  shown  to  provide  a  convenient  tool  for  exposing  and  exploiting  structure  inherent  in  many  model  based 
methods.  It  is  also  shown  that  there  exist  state  space  methods  which  are  robust  to  noise  in  data,  and  to  numerical 
errors.  From  a  computational  point  of  view,  the  methods  are  often  less  complex  than  existing  competing  methods. 
Futhermore  they  only  involve  matrix  operations  which  are  suitable  for  systolic/wavefront  implementation. 


I.  INTRODUCTION 

In  this  paper,  state  space  model  based  methods  for  model  based  signal  processing  are  reviewed  [1-6].  Rational 
modeling  of  signals  followed  by  estimadon  of  the  parameters  of  the  model  has  been  extensively  discussed  in  the  sig¬ 
nal  processing  literature.  Most  often  in  this  context,  a  direct  form  or  transfer  function  parameterization  of  the  models 
has  been  used.  In  this  paper,  state  space  representation  of  these  models  along  with  methods  based  on  these  respresen- 
tations  are  discussed.  The  key  features  of  a  state  space  based  approach  are  summarized  below. 

-1.  The  state  space  framework  provides  a  convenient  mechanism  for  exposing  the  structure  that  is  present  in  model 
based  methods. 

2.  Robust  methods  to  estimate  the  state  space  parameters  of  these  models  are  available.  They  result  in  model  based 
-signal  processing  methods  that  are  numerically  robust,  and  highly  accurate  even  when  noisy  finite  data  records  are 

used. 

3.  The  computational  complexity  of  these  methods  is  usually  less  or  comparable  to  exisung  competitive  methods. 
-Morover,  the  operations  involved  are  matrix  oriented  making  them  suitable  for  systolic/wavefront  implementation. 


H.  STATE  SPACE  MODELING 


In  signal  processing,  digital  filters  or  linear  time  invariant  (LTI)  systems  with  rational  transfer  functions  have 
played  a  central  role.  Most  often  they  are  described  using  constant  coefficient  difference  equations,  i.e. 


y  (n )  =  £  y  (n  - i )  +  £  biu{n-i) 

1=1  1=0 


(1) 


where  u  (n )  is  the  input  and  y  (n )  the  output.  In  the  transform  domain  the  system  is  described  by  its  transfer  func¬ 
tion  H  (z  )  where 


Hiz) 


A(z) 

B{z)’ 
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